Time evolution of local observables after quenching to an integrable model 
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We consider quantum quenches in integrable models. We argue that the behaviour of local 
observables at late times after the quench is given by their expectation values with respect to a single 
representative Hamiltonian eigenstate. This can be viewed as a generalization of the eigenstate 
thermalization hypothesis to quantum integrable models. We present a method for constructing 
this representative state by means of a generalized Thermodynamic Bethe Ansatz (GTBA). Going 
further, we introduce a framework for calculating the time dependence of local observables as they 
evolve towards their stationary values. As an explicit example we consider quantum quenches in the 
transverse-field Ising chain and show that previously derived results are recovered efficiently within 
our framework. 
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Introduction. Recent years have witnessed dramatic 
progress in the study of isolated quantum systems out of 
equilibrium, in particular in systems of optically trapped 
ultra-cold atomic gases. Key to these advances is the 
weak coupling to the environment, which allows the real- 
ization of essentially unitary time evolution on long time 
scales [1-6]. The experimental results have stimulated 
intense theoretical efforts aimed at answering fundamen- 
tal questions such as: Do observables relax to time- 
independent values? What are the principles determin- 
ing these values? How can one describe the relaxation 
towards stationary behaviour? 

There is compelling evidence that nonequilibrium time 
evolution is strongly affected by dimensionality and the 
presence of conservation laws. The experiments of [2] on 
trapped 87 Rb atoms established that three-dimensional 
condensates rapidly relax to a stationary state character- 
ized by an effective temperature, whereas constraining 
the motion of atoms to one dimension greatly reduces 
the relaxation rate of the momentum distribution func- 
tion. These results spurred a flurry of theoretical activity 
aimed at shedding light on the precise effects of inte- 
grability on the nonequilibrium dynamics of many-body 
quantum systems (see [7-34] and references therein). 

So far two basic paradigms have emerged in transla- 
tionally invariant models: at late times subsystems either 
thcrmalize, i.e. are characterized by a Gibbs distribution 
with an effective temperature, or they are described by 
a generalized Gibbs ensemble (GGE) [8]. When the time 
evolution occurs under the action of an integrable Hamil- 
tonian, the GGE is applicable. Questions regarding the 
approach towards the steady state long after the quench 
remain difficult to tackle. Short and intermediate times 
can be efficiently studied by algorithms based on matrix- 
product states [5, 13-15], while numerical methods based 
on integrability have allowed to access arbitrary times in 
finite systems [16-18]. The only cases which have been 



largely understood are noninteracting theories such as 
the transverse field Ising chain [19-24]. 

It is our purpose here to develop an efficient frame- 
work for the description of the out-of-equilibrium dynam- 
ics of a system evolving under an integrable Hamiltonian 
H(h), where h is a system parameter such as an interac- 
tion strength or a magnetic field. Our approach applies 
equally to quantum spin chains and to continuum theo- 
ries like the Lieb-Liniger model. The situation we have in 
mind is that of a quantum quench: a given system is pre- 
pared in the ground state | "J) of the short-ranged Hamil- 
tonian H(ho), which itself may not be integrable. At time 
t = the system parameter is suddenly changed from ho 
to h, and the system evolves unitarily under H(h) for all 
t > 0, i.e. |*(i)) = e _iff W*|*}. Our main focus is the 
calculation of the expectation values of generic, local (in 
space) operators O 
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Examples of O would be products of spin operators lo- 
cated in a finite segment of a spin chain, or density or 
field operators in quantum gases. 

Our main result is to show that in the thermodynamic 
limit L — > oo, at fixed particle density N/L and for a wide 
class of local observables, the expectation value (1) can 
be expressed in a simple way in terms of projections onto 
a single judiciously-chosen representative 'saddle point' 
eigenstate of H{h): 

(*|0(t)|*. 



N 



lim (0(t)) 



lim 

N->oo 



(2) 



In the stationary state, i.e. the limit t — ¥ oo, the averages 
of local observables are simply given by the expectation 
value in the eigenstate 
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We stress that no time averaging is involved in (3), which 
can be thought of as a generalization of the eigenstatc 
thermalization hypothesis [35] to local observables in in- 
tegrable models. We wish to emphasize that the rep- 
resentation (2) offers a dramatic reduction in computa- 
tional complexity as compared to earlier approaches. 

GTBA approach to nonequilibrium evolution. Let us 
consider a quantum integrable model with Hamiltonian 
H solvable by Bethe Ansatz [36]. Let {|4>)} be a complete 
set of eigenstates, i.e. H\&) = The time evolution 

of an arbitrary initial state is then given by 



(4) 



where are constant, complex- valued overlaps 

£i = -ln($|*>. (5) 

Substituting (4) into the numerator of (1) gives a spectral 
representation of the form 

(*|0(i)|*) = e~ £ *- £ *' e^-^'^mOl®'}. (6) 

This double sum over a full Hilbert space basis is a serious 
bottleneck. To proceed further, we look to the thermody- 
namic limit. In the Yang- Yang approach to equilibrium 
thermodynamics [37], a summation over states is recast 
M S ir l functional integral over root densities p 



If) 



\p), £(•••)-►/ 



V\p\e Sp (•••)■ 



(7) 



Here S p is the entropy of all states characterized by a 
given root density and (...) represents quantities with 
well-defined thermodynamic limits. Using (7) once, we 
can formally recast (6) in the form 
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The reason for using (7) only once is that we are inter- 
ested in local operators O. These have the property that 
^ only if both |$) and [$') scale to the same 
distribution p up to microscopic differences [43]. In the 
thermodynamic limit the denominator in (1) becomes 



= J V[p]e- 2Kc{£ » )+s ». 
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and can be evaluated by the method of steepest descent. 
The right-hand side of (9) can be viewed as the partition 
function of an integrable model with "generalized free 
energy" 



T p = 2Re(f p ) - S f 
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Here S p is the usual Yang- Yang entropy of the integrable 
Hamiltonian H{h). In the simplest scalar case, realized 



e.g. in the Lieb-Linigcr model, it takes the form S p = 
N JdX [(p + p h )ln(p + p h ) - p]np- /in/]. The hole 
density p h is related to the particle density p by the ther- 
modynamic form of the Bethe equations 
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where K(X) is a known function for a given integrable 
model. The first term in (10) plays the role of an effec- 
tive energy per temperature and hence acts as the "driv- 
ing term" in a generalized Thermodynamic Bethe Ansatz 
(for details, see [18, 38]). Since the effective overlaps (5) 
are strictly bounded from below, there exists a saddle- 
point at p s , i.e. -j p £ -\p s = [44]. In the thermodynamic 
limit, fluctuations around the saddle point are negligible 
and thermodynamic averages can be calculated with re- 
spect to the energy eigenstatc characterized by p s . Given 
that the expectation values of all local integrals of motion 
in this state are by construction the same as those of the 
generalized Gibbs ensemble corresponding to H(h) and 
\ty(t = 0)), the saddle-point average of local observables 
precisely reproduces the GGE average in the sense of [21] . 
The functional integrals in (8) can be evaluated analo- 
gously: Given that is non-zero only for states 
($| such that w$ — u! p and + £ p are intensive, the 
first term in (8) is dominated by the same saddle point 
p s . The second term is treated analogously. Putting ev- 
erything together we obtain the thermodynamic limit of 
(2). In practice we consider the theory in a large, fi- 
nite volume L (at fixed density N/L) and a particular, 
representative eigenstate |4> s ) that reduces to \p s ) in the 
thermodynamic limit. The corresponding spectral repre- 
sentation is then 
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The gain in efficiency in (12) as compared to the "bare" 
spectral representation (6) is manifold: first of all, a sin- 
gle sum remains; moreover, this sum needs in practice to 
be carried out only over the subset of states with non- 
negligible matrix elements. We note that in our line of 
arguments we have not had to require t to be large. We 
conjecture that (12) describes the time evolution of lo- 
cal observables in the thermodynamic limit, at arbitrary 
times after the quench. Importantly, in the limit t — > oo 
the integrals over the summation over <1> can be carried 
out by a stationary phase approximation. This shows 
that in the stationary state only the expectation value in 
|$ s ) survives in (12), i.e. 



lim (0(t)) 



lim 

N—tOO 



(13) 



The physical content of (12) is as follows. In the thermo- 
dynamic limit (0(t)) is fully determined by states within 
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a basin defined by the saddle point density p s . The latter 
defines particular simultaneous eigenstatcs |$ s ) of all lo- 
cal conservation laws (including the Hamiltonian H{h)). 
At late, finite times (0(t)) is determined by quantum in- 
terference effects due to "excitations" over the state 

An explicit example: the Transverse Field Ising Chain. 
The Hamiltonian of the TFIC is given by 



H(h) 
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where tr" are Pauli matrices at site j of a one dimensional 
chain and we consider J,h>0. At zero temperature and 
in the thermodynamic limit, the TFIC exhibits ferromag- 
netic long-range order along the x-direction for h < 1, 
while it is in a paramagnetic phase for h > 1 [39]. The 
two phases are separated by a quantum critical point in 
the Ising universality class. It is well known that H{h) 
can be diagonalized by combined Jordan- Wigner and Bo- 
goliubov transformations [39] 
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where the single-particle energy is given by £h(k) — 
2Jy/l + h 2 — 2hcosk. Our quench protocol is as follows: 
we prepare the system in the ground state \^) for an ini- 
tial value ho of the transverse magnetic field. At time 
t = we instantaneously change the field from ho to h. 
The state of the system at times t > is obtained by 
evolving with respect to the new Hamiltonian H(h), 



|*(t)) = e -' lH{h)t \^) 



(16) 



The reduced density matrix of a subsystem A at time 
t after the quench is given by PA{t) = TT^p(t) = 
Tr A \^i(t)){9(t)\, in which A is the complement of A. 
For quenches originating in the paramagnetic phase, i.e. 
ho > 1 , the lii symmetry of rotations by 7r around the z- 
axis remains unbroken and it is possible to express PA(t) 
in the form [40] 



PA 
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Here the expectation value is with respect to the state 
|\&(t)) and am and a<i n -\ are Majorana fermion opera- 
tors fulfilling anticommutation relations {aj,a^} = 26j t k, 
which are related to the lattice spins by a Jordan- Wigner 
transformation 
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The matrix W is given by tanh W- = T [41], where 



T jk = Tr [p(t)a k aj] 



kj- 
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In the thermodynamic limit, the correlation matrix is 
given by [22] r 2 „_i,2j-i = r 2 j,2n = /j-n, r 2n -i,2j = 



with 
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where cos A* 



* - 4j2(1 ^ ( : o W )C ° Sfc) - The reduced den- 
sity matrix (17) is Gaussian, and hence multi-point cor- 
relation functions are obtained by Wick's theorem. Con- 
comitantly all local correlation functions in the station- 
ary state are fully specified by the two-point averages 
(19), (20) in the limit t — > oo. So far we have consid- 
ered only the case ho > 1. For quenches originating in 
the ferromagnetic phase, i.e. ho < 1, the reduced density 
matrix pA(t) is not Gaussian [40]. However, as shown in 
[22], /^(oo) is again given by the t — > oo limit of (17). 

Stationary Behaviour. We will now show how to re- 
cover these results in the GTBA framework. The sim- 
plest way to obtain the solution of the GTBA equations 
for the TFIC is to note that the mode occupation num- 
bers constitute conserved quantities [ct k cik,H(h)] = 0. 
Hence the root density in the stationary state is simply 
given by 



p(k) 



(*l4"kl*> 



1 — cos Ai 



(21) 
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and the particle density is D = dk p(k). The corre- 
sponding Hamiltonian cigenstate at density D = 2N/L 
in a large, finite volume is then 



A" 



l$.s 



n<« + - K .|0;/i) 

3=1 



(22) 



where \0;h) is the fermionic vacuum state and the mo- 
menta Kj > are distributed according to the density 
(21), i.e. Kj+i = Kj + l/[Lp(Kj)\. The density matrix 
corresponding to the state (22) is p s = |$ s )($ s | and 
by virtue of the product form (22) it is Gaussian. This 
means that it can be represented in the form (17) and 
is completely determined by its correlation matrix (19). 
The only non- vanishing matrix elements are 



1 2l-X,2l-2n 



i x ^ e 

k 



nk (h-e lk )(l-2S k , Kj ) 
\/l + h 2 - 2h cos k 



(23) 

Turning the sum over momenta into an integral by means 
of the Euler-Maclaurin sum formula, we find that T s = 
r(oo) and hence lim^oo pA{t) = Ps,A- This proves that 
the GTBA formalism reproduces the correct stationary 
state for the reduced density matrix for any finite subsys- 
tem in the thermodynamic limit, and hence for all local 
correlation functions. 
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Relaxation behaviour. Our general formalism suggests 
that the time evolution of local (in space) operators is 
given by (2), where the state |<& s ) is defined in the previ- 
ous paragraph. We now demonstrate the validity of (2) 
for any local operator O in the case where the quench 
originates in the paramagnetic phase, such that the Z2 
symmetry is unbroken. The proof is as follows: we start 
by defining two density matrices p{t) = \^f(t))(^(t)\ and 
p s (t) = \$ s (t)){&s(t)\- Crucially, both of these density 
matrices are Gaussian as a Wick's theorem holds for av- 
erages calculated with respect to both of them. The right 
hand side of (2) can be written in the form 



Tr [p s (t)p(t)Q} 
2Tr[ Ps (t)p(t)} 



+ P^ Ps 



= Tr \p(t)0] 



(24) 



Because each term in p is a product of two Gaussian 
density matrices, it is Gaussian itself, and hence fully 
characterized by its correlation matrix 



T jk = Tr [p(t)a k aj] - 8j, k . 



(25) 



The two-point functions in (25) are easily calculated, 
and as shown in the supplementary material we have 
lim7v^ oo r(t) = r(i). This proves that for any local 
observable O (such as products of spin operators con- 
tained in a finite subsystem) in the thermodynamic limit 
liniA^oo Tr [p(t)0] = Tr [p(t)0], and establishes eqn (2). 

Our line of arguments breaks down for quenches orig- 
inating in the ferromagnetic phase, because the den- 
sity matrix p(t) is no longer Gaussian [40]. In order 
to verify the validity of (2) in this case, we have an- 
alyzed the relaxation of the order parameter one-point 
function (^f (t)\a^ (t)) for quenches with ho,h < 1 in 
the regime where the density of excitations of the post- 
quench Hamiltonian H(h) in the initial state is small, i.e. 
{^\a\a k \^) < 1. The result for Jt > 1 in this case is 
[211 



!*(*)>« (l-ft 2 )« exp 



-2t f d ±e' h {k)K\k) 



(26) 

This result is recovered from (2) in a very efficient way 
as follows. Taking into account boundary conditions in a 
large, finite volume, the appropriate form of (2) for the 
order parameter expectation value is 



(tt(t)|of|*(t))=Re 



NS 



s/NS 



(27) 



Here R/NS correspond to pcriodic/antiperiodic bound- 
ary conditions on the fermions respectively, see e.g. [21]. 
As shown in [21], the state r(\E'(£)| can be written as a 
linear superposition of energy eigenstates with n pairs of 

fermions R (*(t)| = £n=o *<**(*)!• The late-time be- 
haviour of (27) is determined by states with N pairs, i.e. 
the term with n = N. Retaining only this contribution, 



and using the known form of matrix elements of the or- 
der parameter [42], one readily obtains the result (26) by 
means of the techniques developed in [21]. More details 
are provided in the supplementary material. 

Conclusions. We have argued that averages of local 
observables in the steady state reached long after a quan- 
tum quench to an integrablc model can be described as 
expectation values with respect to a single cigenstate of 
the Hamiltonian (eqns (3) and (13)). This cigenstate 
can be constructed by means of a generalized Thermo- 
dynamic Bcthe Ansatz and corresponds to a saddle point 
of an appropriately defined "free energy" . Going further, 
we have shown that the time evolution of local observ- 
ables is governed by states in the vicinity of this saddle- 
point through eqns (2) and (12). The spectral repre- 
sentation (12) allows to identify the physical mechanism 
underlying the relaxation for a given observable at late 
times. Our approach paves the way for analyzing quan- 
tum quenches in interacting integrable models and appli- 
cations to the Lieb-Liniger and the sine-Gordon model 
are in progress. An interesting question concerns the 
range of initial states that can be analyzed within our 
framework. In quenches of the kind we have focussed 
on here, the initial state is translationally invariant and 
characterized by low entanglement. Whether these are 
necessary requirements and whether the regularity as- 
sumptions underlying our GTBA analysis hold for any 
such states are important points for further investigation. 
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Supplementary Material 

Calculation of the correlation matrix T (25) 

For quenches originating in the disordered phase the 
Z2 symmetry is unbroken and no subtleties emerge with 
regards to boundary conditions in a large, finite volume 
L. In particular all momenta in the formulas below can 
be taken in the NS sector, i.e. p n = ^-(n + 1/2), where 
—L/2 < n < L/2. The state ^(i)) can be expressed in 
terms of eigenstates of H(h) as [21] 
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\0;h), (28) 



where K{k) = tan(Afe/2), and a\ obey canonical 
anticommutation relations, and ak\0',h) = 0. The time 
evolved state |$ s ) (22) is 

N 

='[[e- 2is ^ t ai j al Kj \0;h) . (29) 
j'=i 

Our goal is to determine all two-point matrix elements of 
the form (t)\a,jak\& s {t)) , where aj are the Majorana 
fermion operators defined in (18). The mode expansion 
of the a,j is 
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A straightforward calculation using the mode expansion 
and a p \0; h) = gives 
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In the limit TV , L — > oo with D = 2N/L fixed we can 
turn sums into integrals using the Euler-Maclaurin sum 
formula 



i N r 

7 E/(«i) = dkp{k)f{k) + 0{L- 1 ), (33) 

where p(fc) is the root density (21) characterizing the 
distribution of the momenta {nj}. Using that K(k) = 
tan(Afc/2) we find 

JToo <¥^) = 5 » -° + ' (34) 

where /; is given in (20). By taking the complex conju- 
gate of (34), we find 

A^oo = Sn -° (35) 

Given that /„ is purely imaginary, we conclude that 

r2n+i,i = Tr [p(t)aia 2n +i] - <5«,o = /-«• (36) 

By analogous calculations one finds that 

(*(t)\a2a2n+2\<S>s(t)) 
hm = d„. + Jn, 



(*(t)|a 2 a 2n+1 |$ s (t)) 
hm . . ; = q n . 



(37) 



Taking the complex conjugates of eqns (37) and using 
that both /„ and g n are purely imainary, we obtain 

T2n+2,2 = Tr [p{t)a2Cl2n+2] — 5n,0 = fn , 

T2n+i,2 = Tr [p(t)a 2 a 2n +i] = g n - (38) 



Finally, using translational invariance we conclude that 



r — r 

L 13 ~ 1 n i 



(39) 



where Tij is given in (19), (20). 



Order Parameter One-Point Function 

As we are interested in the order parameter, we need 
to carefully keep track of boundary conditions. We follow 
the conventions summarized in the appendix of Rcf. [21]. 
The state after the quench is given by [21] 



!*(*)> 



|tt(f))R + |tt(*))NS 

V2 



— r l + ie-^^Klk)^ 

|*(*)>. = -*■ fc |0;fe) a ,(40) 

o<ke. ^{l + KHk)) 
where a = R, NS. The appropriate form of (2) reads 



(*(t)|o?|*(t)) 



lim Re 

JV->oo 



(*|of(t)|*,)NS 



s/NS 



(41) 



The state R,(^(t)| can be written as a linear superposition 
of energy eigenstatcs [21] with n pairs of fermions 



L/2 

R (*(t)\ = Y^ K (* n (t)\ 

n=0 



(42) 



Inspection of the resulting matrix elements suggests that 
at late times the dominant contribution to (41) arises 
from terms with TV pairs. The corresponding contribu- 
tion is 



(y N (t)\af |$ fl (t)) NS _ (-1 



NS 



NS 



TV! 



^2 k(Pi,-Pi,---Pn,-PnM\ki,-ki, 



0<pi 



,PN 



\ TJ K(Pj) 2it(e v .-e K .) 

(43) 



where e p = £h(p)- The matrix element of erf in (43) is Here £ = |1 — h 2 ] 1 / 4 , e PA = (e p + e q )/2 and c VA = 
given by [21, 42] [cosp — cosq]" 1 . As a function of pj it exhibits either 

simple or double poles when pj — > k m , and the lead- 
(— l) Ar 2 2Ar y / ? TT 4 2 tt shiKj smpj ing contributions at late times arise from the vicinities 

T2N ll e t 3 ,Pl C K 3 ,Pl 11 F 2 F 2 

3,1 3 = 1 «i * 

X II £ 4 £ 4 \2 — • ( 44 ) 

j<j> Kj, K j< Pi >P}> k j Pj -Pj' 
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of these singularities [21]. We proceed by rewriting the 
"pole-factor" as the square of a Cauchy determinant 



JV 



i 



3,1=1 j < ji C Vj>Pj' C ^i^j' 



Ci 



COS Pi — COS Kj 



(45) 



By expressing the determinant as a sum over permuta- 
tions using the standard definition, we can recast V in 
the form V = ^jlo^ ^N-ij, where Vk has k double poles 
when viewed as a function of the cosp r . The most sin- 



gular terms are (x r = cosp r , y s = cosfc s ) 



Vn - E II t x .\ p y 



N-2 



9 E E II ( Xj - yPj y 



P^Sn r^s j^r,s 

l 



(x r - yp r )(x r - yp s )(x s - yp r ){x s - y Ps ) 



(46) 



The leading contribution to (43) at late times is obtained 
by retaining only the most singular part Vn 



(\Mt)|<7?|<Mi)) NS „ 2 2W ^ 



N 



NS 



NS 



e 



0<pi 



K(j>j) sin re j sinpj 



n 

3,1 



n 

i<i' 



(47) 



p 3 -p,- 



We can now carry out all sums over pj using (D.l) of 
Ref. [21]. Retaining only the leading contributions at 
late times (and asymptotically large L) we find 



R(*Jv(*)l^l* fl (*)>NS 



N 



NS 



s/NS 



3=1 

V^exp 



4f 
7 £ 



AJ. N 

3 = 1 



(48) 



Turning the sum over the Kj into an integral this becomes 

R<*w(*)|of|*.(*)>NS ^ ^-t/o-^ (l-cos A,) _ (49) 



NS 



Observing that 



1 — cos At 



2K 2 (k) 
1 + K 2 (k)' 



(50) 



and noting that the asymptotic result (49) is purely real, 
we recover the result obtained in [21] by summing the 
leading late-time divergences in the form factor expan- 
sion for the 1-point function to all orders. We stress that 
the current calculation is significantly simpler. Sublead- 
ing contributions at late times can be obtained by taking 
into account V^-ij with j > 1. These will be reported 
elsewhere. 



s NS 



